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A novel method for simulating the statistical mechanics of molecular systems in which both nu- 
clear and electronic degrees of freedom are treated quantum mechanically is presented. The scheme 
combines a path integral description of the nuclear variables with a first-principles adiabatic de- 
scription of the electronic structure. The electronic problem is solved for the ground state within a 
density functional approach, with the electronic orbitals expanded in a localized (Gaussian) basis 
set. The discretized path integral is computed by a Metropolis Monte Carlo sampling technique 
on the normal modes of the isomorphic ring-polymer. An effective short-time action correct to 
order r 4 is used. The validity and performance of the method are tested by studying two small 
Lithium clusters, namely Li4 and Lig~. Structural and electronic properties computed within this 
fully quantum-mechanical scheme are presented and compared to those obtained within the classical 
nuclei approximation. Quantum derealization effects turn out to be significant as shown by the fact 
that quantum simulation results at 50 K approximately correspond to those of classical simulations 
carried out at 150 K. The scaling factor depends, however, on the specific physical property, thus 
evidencing the different character of quantum and thermal correlations. Tunneling turns out to be 
irrelevant in the temperature range investigated (50 K to 200 K). 
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I. INTRODUCTION 

Light atoms, such as H, He, Li or Be, cannot very of- 
ten be treated as classical particles, particularly at low 
temperatures. As temperature decreases and the ther- 
mal de Broglie wavelength increases, the quantum char- 
acter emerges, and a description in terms of classical co- 
ordinates and momenta breaks down. The most obvious 
manifestation of the quantum character of light atoms is 
a large zero-point-energy (ZPE). A particle of mass m 
in a harmonic potential with characteristic frequency ui 
will have a ZPE of hw/2 and an associated spatial der- 
ealization of Ax = yjh/muj. For instance, a proton in 
a typical bonding environment such as a H-0 bond or 
a H 2 molecule, will have a ZPE of 0.15 to 0.25 eV and 
Ax between 0.2 and 0.3 rA. This represents a sizeable 
effect which could be decisive in stabilizing a particular 
crystalline structure for a solid, or the ground state con- 
figuration of a molecule or a cluster. An even more inter- 
esting manifestation of quantum effects is the possibility 
that these light nuclei can tunnel across potential energy 
barriers, thus exploring classically forbidden regions of 
configuration space and giving rise to a variety of inter- 
esting quantum effects such as temperature-independent 
diffusion, exotic ground states, resonances in ion-surface 



scattering, and fluxional molecules. Signatures of quan- 
tum effects can also be observed in low-energy atomic 
collisions, or in proton-transfer reactions in the gas and 
condensed phases. 

To date, most studies that consider the quantum char- 
acter of atomic nuclei are based on an empirical descrip- 
tion of the interatomic interactions, or otherwise con- 
sist of extending and/or correcting a posteriori the re- 
sults obtained within a classical nuclei approximation. 
Classical potentials are frequently not transferable from 
one environment to another, and are ill-suited to mod- 
eling phenomena involving significant electronic density 
redistribution, as in the making and breaking of chemical 
bonds. The natural route to overcome these limitations 
is to describe the interactions at a first-principles level, 
i.e. by including explicitly the electronic component in 
the description of the system. The recent development of 
such schemes, which address the question of the interplay 
between electronic structure and the quantum nature of 
light nuclei, has only very recently began to be realized, 
thus opening a fascinating field with important implica- 
tions in many branches of physics, chemistry and biology. 

Since small clusters usually exhibit rich ..landscapes of 
isomeric forms within narrow energy bandsEJ, they consti- 
tute good systems for studying the effects of quantum de- 
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localization and associated tunneling behavior. Lithium 
clusters are particularly interesting because in addition 
to the small atomic mass, they are bound by metallic 
many-body interactions which cannot be adequately pep- 
resented by means of classical interatomic potentialscl. 

In the remainder of this Introduction we outline the 
methodology that we have developed to study this class 
of problems, and review the present understanding of Li 
clusters which is the test system for our method. In Sec- 
tion II we introduce the theoretical framework of our ab 
initio path integral approach and discuss the approxima- 
tions involved. Section III is devoted to the details of 
the practical implementation of the path integral Monte 
Carlo (PIMC) and the electronic structure methods. In 
Section IV we present the validation of the electronic 
structure calculations, zero-temperature geometries and 
electronic properties of Li4 and Li^~ clusters. The re- 
sults of our simulations for the classical and quantum 
Li4 and Lijj~ clusters at finite temperatures are presented 
in Section V. Section VI contains our conclusions and an 
assessment of the potential of this novel simulation tool. 



A. Methodological aspects 

The goal of the present work is to introduce a novel 
computational technique for studying the statistical me- 
chanics of isolated systems like clusters and molecules 
containing light atoms. Our approach combines an 
imaginary-time path integral description of the nuclear 
degrees of freedomB with a first-principles density func.- 
tional (DFT) description of the electronic structurea. 
Since the natural choice for investigating isolated systems 
is to use a localized basis set fpf the electronic orbitals, 
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we adopt a Gaussian basis setcl. In the present imple- 
mentation the electronic structure is computed at the 
all-electron level, i.e. explicitly including core electrons. 
The sampling of the path integral is implemented using 
Monte Carlo (MC) techniques. The electronic energy is 
minimized for each nuclear configuration, and the MC re- 
jection step is also performed using the energy calculated 
at the same level of sophistication. 

Other schemes along these lines have been recently 
proposed by Marx and ParrinclloO, and Cheng et al.Q. 
At variance with our approach, these two methods use 
Molecular Dynamics (MD) for sampling the path inte- 
gral, a plane-wave (PW) expansion for the electronic or- 
bitals and treat the electron-ion interaction at the pseu- 
dopotential level. PW expansions with periodic bound- 
ary conditions are more appropriate for extended systems 
such as solids or liquids though they can be adequately 
modified to deal with isolated systemsQ. 

The evaluation of real-time path integrals, which would 
include the full dynamical information, is an extremely 
difficult numerical task because the integrand is a rapid 
oscillatory function of the path. In the imaginary-time 
framework the statistical weight is real and positive def- 



inite, so that the integrals are well-conditioned, but the 
price is that the dynamics is not directly accessible. In 
the absence of real-time dynamical information, it is 
not particularly advantageous to sample the integral us- 
ing MD in place of MC techniques. In particular, the 
MD technique requires elaborate thermostatting mecha- 
nismsu in order to overcome ergodicity problems in the 
sampling of the quasi-harmonic degrees of freedom that 
appear in the path integral formulation (see Section II). 
In the present method we propose a Metropolis MC sam- 
pling technique on the normal modes of the polymer, 
which has less severe ergodicity problems, and is more 
convenient from the point of view pf efficient evaluation 
of the path action (see Section 111)0. Moreover, the MC 
strategy is easy to interface with any electronic structure 
code. 



B. Small Lithium clusters 

Structural and electronic properties of Li„ and Li+ 
clusters have been systematically investigated for n — 
2 — 9 by paeans of ab initio configuration-interaction cal- 
culationsEl A key observation was the existence of sev- 
eral isomers of comparable energy. Recently, ab initio 
classical MD simulations at the Hartree-Fock (HF) and 
non-local density functional (NLDF) level were carried 
out in order to explore the cluster dynamics as a function 
of temperature, to analyze isamexization reactions, and 
to study the melting transitionOilj. The extent to which 
different levels of first-principles calculations are reliable 
for describing small Li clusters has been very recently 
discussed by Rousseau and MarxEj, who concluded that 
either ab initio MP2 calculations, or gradient corrected 
NLDF provide a reasonable potential energy landscape, 
while HF and LDF are inadequate. 

This aspect is very important because the energy land- 
scape determines the probability with which the cluster 
visits different possible geometries. In general, the abil- 
ity to jump from one minimum to another will depend 
on the extent of both thermal and quantum fluctuations. 
Fig. 1 shows three typical situations: At low temper- 
atures quantum derealization is the dominant mecha- 
nism. If the ZPE is higher than the barrier between the 
two minima (a), or smaller than the barrier but larger 
than the energy difference between them (b) (tunneling 
regime) , the ground state wave function will sample both 
configurational minima. If, instead, it is smaller than 
the energy difference between the two minima (c), then 
the ground state will be basically unchanged from the 
classical one, though all structural distributions will be 
broadened by quantum derealization. Thermal excita- 
tions can promote a classical-like situation like (c) to a 
thermally-assisted tunneling regime. 

The above picture is valid as long as the electronic 
ground state is non-degenerate, and this is realized for 
clusters with closed electronic shells. Open-shell clus- 
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ters undergo Jahn- Teller distortions, and are character- 
ized by degenerate ground states and pseudo-rotations. 
An adequate treatment of this situation would require 
the introduction of the-, concepts of conical intersections 
and geometric phaseal3. Though these phenomena ate 
of greatinterest, and have indeed been detected in LiatZl 
and Lislia, we have preferred to avoid this extra compli- 
cation at this stage in the development of our method. 

Quantum nuclear effects in Ii_clusters were first ad- 
dressed by Ballone and MilanO, who studied magic- 
number clusters (20, 40 and 92 atoms) by means of PIMC 
simulations using a simple, jellium model description of 
the electronic component. Their most interesting obser- 
vation was the existence of a large number of isomers 
differing only in the location of the outermost atoms, 
such that quantum tunneling amongst these isomers led 
to a fluid-like behavior. A very recent ab initio path inte- 
gral MD study of Lis and L120 clusters by Rousseau and 
Marxo shows that this picture does not hold when the 
description of the electronic component is improved, thus 
demonstrating the necessity of going beyond the jellium- 
type description. 

For our study we have chosen the following two exam- 
ples: L14 has a single isomer and is well described by a 
quasi-classical picture. Li^", has two isomers at an en- 
ergy difference of about 150 K, and two other isomers 
at higher energies (see section 5). The zero-point-energy 
is of the order of 100 K, so that LI5 could constitute a 
good candidate for exhibiting significant quantum effects 
in terms of the sampling of configurational minima. The 
results presented in section 5 will show that, in spite of 
this, Lit actually behaves in a very similar way to L14. 



II. THE AB INITIO PATH INTEGRAL 
PARTITION FUNCTION 

The statistical mechanics of quantum many-body sys- 
tems can be formulated in terms of the two-point density 
matrix, or imaginary-time propagator: 



P (K, R, p) =< R| exp(-(3H)\R' >, 



(1) 



whose trace is the partition function Z. (3 = 1 /ksT is the 
inverse temperature. The path integral representation of 
the density matrix is given by 



p(R, R',/3) = J T>[R{u)} e" 



S[R(u)] 



(2) 



where R(it) represents the configuration of an ./V-body 
system as a function of imaginary time u. The range of u 
is from to (3Ti and the paths considered are restricted to 
those beginning at R(0) = R and ending at H((3fi) = R'. 
The partition function can be similarly expressed as a 
path integral with contributions from all possible cyclic 
paths for which R(0) = H(/3h). T)[R(u)) represents the 
differential element for all paths. The Euclidean action, 
5[R(ti)], associated with a path is defined as 



S[K(u)} = - 



■in 
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V ~2 



rfR 



du 



V 



R(u) )du. (3) 



The first term corresponds to the kinetic energy con- 
tribution to the action, with to the mass of the parti- 
cles. The generalization to heterogeneous systems, i.e. 
composed of two or more species of different masses, is 
straightforward. 

In order to devise a feasible computational scheme, the 
path integral is typically discretized by representing the 
cyclic paths as a finite set of 37V-dimensional configu- 
rations, R^, at equispaced points in imaginary time be- 
tween and f3h. The degree of discretisation is referred 
to as the Trotter number, P. The short-time or high tem- 
perature propagator, p(Ri, Ri+i; can be evaluated 
semi-classically at different levels of approximation. The 
contribution of the kinetic energy term to the short-time 
action is written in terms of a first-order finite difference 
between configurations on adjacent time slices, while the 
short-time integral of the potential energy together with 
higher-order corrections to the kinetic energy, is replaced 
with an effective, quantum-corrected potential V e ff (R)- 
The sum of the two terms is referred to as the effective 
action. Therefore, the expression for the partition func- 
tion of N interacting, distinguishable quantum particles 
with Trotter number P is given by: 



-■NP 



mP 
2ttH 2 P 



3NP/2 




exp 



\ ^ i—l i—l / 



(4) 



According to the level of approximation of the effective 
action, the number of slices P needed to achieve conver- 
gence in the partition function can be small enough that 
the problem is tractable, or large enough that the evalua- 
tion of the multidimensional integral becomes a hopeless 
task. It is therefore important to use the best possible 
effective action compatible with the computational com- 
plexity involved in its calculation. The simplest one, or 
primitive approximation, replaces the effective potential 
by the bare potential, which is equivalent to an end-point 
approximation for the short-time integral: 



V 



R(u) 



du 



V(R) + V(R') 



(•5) 



and is correct only to order r 2 . At the other extreme, 
the pair-action approximation provides a very accurate 
technique when the full, many-body potential can be rea- 
sonably approximated with a sum of pair potentials. This 
scheme has been very effectively exploited to investigate 
the properties of liquid and superfiuid He down to tem- 
peratures of about 1 KE3, a task that would not have 
been possible using the primitive action. 

As mentioned in the introduction, classical interaction 
potentials are computationally fast, but very often un- 
reliable. Realistic interaction potentials can instead be 
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obtained from more expensive first- principles techniques. 
In these latter, the electronic degrees of freedom are ex- 
plicitly included in the Hamiltonian description of the 
system: 

W(R,r) = f n + f c + y cc (r) + T/ on (r,R) + K n (R), (6) 

where r and R are the electronic and nuclear coordinates, 
T and V stand for kinetic and potential operators, while 
subscripts e and n indicate electronic and nuclear compo- 
nents, respectively. The path integral representation for 
the partition function could then be developed using thje 
coordinate basis for both the electrons and the nucleiE3. 

However, standard electronic structure calculations are 
carried out in a wave function representation by resort- 
ing to the adiabatic separation of nuclear and electronic 
motion. It is therefore more convenient to expand the 
electronic component in the adiabatic basis set where 
electronic wave functions \<p a > and total energies E a (R) 
arc obtained by diagonalising the electronic Hamiltonian 
T e + V ee (r) + t4n(r,R) + Kn(R). If t = (3/P, then the 
discretized partition function reads: 

- Z F=X)'"X] /■"/ II ( P««,a i+1 (Ri,Ri+l)T) d&i 
ai a P J i=l \ , 

(7) 

where 

p a , 7 (R,R',r) =< R| < <f> a \ exp(-T7i) ^ > |R' > . 

(8) 

Then, in the spirit of expression (Q), the short-time prop- 
agator can be written as: 

p a , 7 (R,R',r) =< R| < 4> a \ exp(-rT n ) |</> 7 > |R' > x 

exp[~[E e Jf(K) + E^(K')]), (9) 

where E e Ji(R) is an effective potential which derives di- 
rectly from the electronic structure. In the primitive ap- 
proximation, which is what-Jias been used in other ab 
initio path integral methodstffl, the effective potential is 
simply the total energy corresponding to adiabatic state 
a of the electronic Hamiltonian, i.e. E a (R). 

The nuclear kinetic energy operator can give rise to 
non-adiabatic coupling matrix elements between adia- 
batic eigenstates. If these are negligible, but more than 
one Born-Oppenheimer (BO) surface is occupied due to 
thermal excitations, then the partition function will split 
into independent manifolds indexed by the BO electronic 
eigenstate. In the absence of degeneracies in the ground 
electronic state, the energy differences between electronic 
eigenstates are typically orders of magnitude larger than 
reasonable thermal kinetic energies. Consequently only 
the ground electronic state contributes to the partition 
function and the electrons only enter at the level of re- 
placing the total potential energy with the ground state 



first-principles effective potential Eq (R) in equation 
(I): 

=p (-|£X>. - R-.«) 2 - P t.^"(Ro) ■ (io) 

This is the partition function that will be evaluated us- 
ing Monte Carlo techniques. The expression for Zmp 
can be interpreted as the partition function of TV classi- 
cal polymers, each of P monomeric units or beads, with 
adjacent beads linked by harmonic springs with force con- 
stant mP/(3h 2 . Beads on two separate cyclic polymers 
are coupled by the interaction potential only if they lie 
on the same time slice. 

Let us remark that the computation of excited elec- 
tronic states is an open issue in density functional the- 
ory, and it is known that the usual approximations to ex- 
change and correlation, like the LDA and even NLDF do 
not provide reliable excitation energies. Excited states 
could be calculated properly using very high-level ab 
initio methods. However, since we have developed a 
methodology in which the electronic degrees of freedom 
are dealt within DFT, we are for the time being not in a 
position to incorporate non-adiabatic couplings and ex- 
cited electronic manifolds. 



III. PRACTICAL IMPLEMENTATION 
A. Path integral Monte Carlo 

1. Effective action 

We use a discretized time representation in which a 
path is described as a set of configurations, {R;},i = 
1, • ■ ■ P, at P-equispaced points in imaginary time. The 
effective short-time propagator for two adjacent points 
along the path has been evaluated to fourth-order ac- 
curacy in t = J3JP, so that the first-principles effective 
potential readstJ: 

*"(*>- *Ca.) + (££0 g(^)) 2 , 

(11) 

where Xjj is the 3-dimensional vector of the coordinates 
of particle j in slice i, such that R^ = (xii,Xi2, ■ • ■ , Xjjv). 
This form of the effective action leads to an error of the 
order of P((3/P) 5 in the partition function and allows 
us to significantly reduce the Trotter number required 
for convergence as compared to the primitive action (see 
Section V). 
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The quantum correction to the potential requires the 



evaluation of the first-principles forces F,- 



-dE /dx.j 



in the ground electronic state. The cost of this oper- 
ation is of the same order of magnitude as the rest of 
the electronic structure calculation, at least in most ab 
initio or density functional schemes. In contrast, second 
and higher-order derivatives are sufficiently expensive to 
evaluate that the computational advantages of a more 
accurate short-time action are lost. 



2. Normal modes sampling 

The above expression for the partition function can 
be directly used to set up a Metropolis MC simulation 
scheme by assigning the appropriate Boltzmann weight to 
each configuration of the N x P-dimensional isomorphic 
classical system. However, as quantum effects increase, 
the degree of discretisation must be increased to main- 
tain accuracy. Since the harmonic force constant between 
adjacent beads on the quantum polymer is mP/(3h 2 , in- 
creasing the Trotter index results in increasingly stiff 
harmonic links and the computational problem of ensur- 
ing the ergodicity of the Metropolis walk becomes in- 
tractable. An intuitively appealing and computationally 
simple way for circumventing this difficulty comes from 
considering the normal modes of the quantum polymeria. 
In the absence of an interaction potential, all Cartesian 
degrees of freedom of the system are decoupled. For a sin- 
gle degree of freedom the harmonic intra-polymer poten- 
tial is given by V p = (mP/2(3h 2 ) Yli=i( x l — xi+i) 2 - Diag- 
onalization of the second derivative matrix of this poten- 
tial leads to the normal coordinates, {Qk}, k= 1, • • • , P, 



Qk = (l/y/P)^2xiexp(2irikl/P). 



(12) 



i=i 



In the normal mode representation, the kinetic energy 
contribution to the path action is 



0h 



m / dx\ 
~2\~du) 



2 9™ p _L_ 

du = —^J2\Qk\ 2 sin 2 (7Tk/P). (13) 



The zero-frequency mode (k = P) corresponds to mo- 
tion of the center-of-mass of the polymer and makes 
no contribution to the kinetic energy. All other nor- 
mal modes would be Gaussian distributed with variance 
o~ 2 = (3h 2 1 AmP sin 2 {irk / P) if they corresponded to free 
particles. The potential energy term couples these nor- 
mal modes and cause distortions from the free-particle 
distribution. The low-frequency modes correspond to 
large, collective motions of all beads of the polymer, while 
the high-frequency modes cause small, local path fluctu- 
ations. The normal modes are then used as Metropolis 
variables and the displacements scaled according to the 
Gaussian dispersions associated with each normal mode. 



3. Observables 

The canonical ensemble average of an observable O is 
jiven by 



(O) = Tr{pO}/Tr{p}, 



(14) 



where O is the corresponding quantum mechanical op- 
erator. If the operator O is diagonal in the coordinate 
representation, then 



(O) = J dR O(R) p(R,R;/3). 



(15) 



Our MC strategy samples configurations R with prob- 
ability proportional to p(R, R;/3), such that equilibrium 
averages can be readily estimated via discrete summa- 
tions. Dynamical variables that can be related to the 
partition function are also straightforward to obtain us- 
ing thermodynamic estimators. With our choice of the 
short-time action, thermodynamic estimators of the to- 
tal energy (E), the kinetic energy (K) and the potential 
energy (V) are given by 



(V) - (U 2 ) + 2(U C ) 



(K) = (U 1 ) + (U c ) 



(E) = (U 1 ) + (U 2 )+3(U C ), 



where 



3NP ' mp(Rj _ Ri+i)2 

i=l 



2(3 



2(3 2 h 2 



1 p 

U2 = p J2 E ^ 



i=l 



c P 24mP 2 2^ 2^ \ ft*- 

i=i j=i y 



(16) 
(17) 
(18) 

(19) 
(20) 
(21) 



B. Electronic structure calculations 

The calculation of the electronic energies is carried 
out within the framework of DFT. For a given nuclear 
configuration, Kohn-Sham single-particle equationstJ are 
solved self-consistently for the electronic density, and the 
total energy and forces are computed accordingly. Kohn- 
Sham orbitals are expanded in a Gaussian basis set. 

The electronic density is also expanded in an additional 
Gaussian basis setc3. The coefficients for the fit of the 
electronic density are computed by minimizing the error 
in the Coulomb repulsion energy. The use of this pro- 
cedure results in an important speedup, since the cost 
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of evaluating matrix elements reduces from 0(N 4 ) to 
0(N 2 M) (TV is the number of functions in the orbital 
set, and M the number of functions in the auxiliary set, 
typically of size comparable to N). 

Matrix elements of the exchange-correlation potential 
are calculated by a numerical integration scheme based 
on grids and quadratures proposed by BeckecZI. During 
the self-consistency cycle, the integration is performed on 
a set of coarse atom-centered, spherical grids. At the end 
of the self-consistent procedure, the exchange-correlation 
energy is evaluated using an augmented, finer grid. This 
strategy of combining coarse and fine grids results in a 
considerable gain in computational efficiency which is 
very important because this part is one of the main bot- 
tlenecks of the calculation. 

The exchange-correlation term is described a gradient 
corrected NLDF level. Correlation is given by the pa- 
rameterization of the homogeneous electron gas of Vosko 
et al.E3 supplemented with the gradient corrections pro- 
posed by PerdewcJ. Gradient corrections to the exchange 
term are taken from BeckeEJ. 

The first derivatives of the energy with respect to the 
nuclear coordinates, required by the fourth-order effec- 
tive action, are evaluated by taking analytical derivatives 
of one-electron and Coulomb terms, while the exchange- 
correlation contribution is obtained by numerical inte- 
grationEl. 



IV. GROUND STATE PROPERTIES OF 
CLASSICAL LI 4 AND LI+ 

A. Validation of the basis set and optimized 
geometries 

We have analyzed five different basis sets for Li4 and 
Li5 + r _.The first one (labeled 1) is the standard 3-21G 
basis£j. The second set (labeled 2) is a double zeta plus 
polarization basis set, optimized for DFT calcultions — . 
The third one (labeled 3) is the standard 6-311G basistS, 
the fourth set (labeled 4) is the 6-311G set augmented 
with a polarization function (6-311G*), and finally the 
fifth set (labeled 5), consists of a large uneontracted ba- 
sis set (13s/9p/ld), proposed by DunningEil. The calcu- 
lations performed with basis sets optimized for standard 
ab-initio calculations (labeled 1, 3, 4, and 5) have been 
carried out using an uncontracted auxiliary basis set with 
a scheme (7s/3p/3d), as proposed in Ref. |3^ and The 
calculations performed with basis set 2 have been car- 
ried out using an auxiliary basis set proposed in Ref. |3^, 
with a scheme (7s/2p/ld). Full geometry optimizations 
without symmetry constraints have been performed in all 
cases. 1 — 1| — 1 

In agreement with previous worklj'Ej, only one stable 
minimum with a rhombus geometry has been found for 
L14 (see Fig. 2), while for Li5 + we found four stable lo- 
cal minima (also shown in Fig. 2). The highest energy 



isomer, i.e. isomer I, consists of two triangles which lie 
on the same plane, joined by a shared central atom. The 
second isomer, i.e. isomer II, is similar to isomer I, but 
now the triangles lie on perpendicular planes. The third 
isomer, i.e. isomer III, has C2 symmetry, and can be de- 
scribed as an isosceles triangle plus a dimer. The dimer 
is located perpendicularly to the plane of the triangle, 
close to its shortest side. The fourth isomer, i.e. isomer 
IV, has a trigonal bipyramidal structure. It should be 
pointed out that neither isomer II nor isomer III have 
been reported in earlier works, probably because of sym- 
metry constraints used during the geometry optimization 
procedure. The basis set dependence of the binding ener- 
gies is shown in Figure 3. Relevant structural parameters 
for all basis sets considered here are given in Table I. Our 
results for basis set 4 agree with those reported for the 
same functional and basis set in Ref. [IB] for Li4 and for 
isomers I and IV of Li?. 

The data in Figure 3 and Table I show that calculations 
performed using basis set 2 yielded results that deviate 
considerably from the ones obtained using the very large 
basis set 5, which can be considered as almost converged. 
Even if the errors in computed bond lengths and binding 
energies are not too large (about 5 % for bond lenghts 
and 5-10 % for binding energies), calculations carried out 
using basis set 2 fail in reproducing the energy sequence 
of isomers of Li5 + . This can be ascribed to the fact that 
the description of the p-shell using this set is rather poor, 
since it contains only one function. On the other hand, 
calculations performed with basis sets 1, 3 and 4 yield the 
same energy ordering for isomers of Li5 + and structural 
results within 2-3% from the almost converged, basis set 
5, values. _ 

Basis set superposition error (BSSE)El calculations for 
Li 4 yield 14.36, 1.00, 0.96, 2.34, and 0.13 kJ/mol, for 
basis sets 1, 2, 3, 4, and 5, respectively. It is clear that in 
order to reduce BSSE, and obtain meaningful interaction 
energies, a better description of the p-shell than the one 
provided by the small basis set 1 is required. 

TABLE I. Selected geometrical parameters of L14 and Lis 4 " 
(in rA). Atoms are labeled as in Figure 2. 



System 




Set 1 


Set 2 


Set 3 


Set 4 


Set 5 


Li 4 




di2 


2.658 


2.678 


2.638 


2.625 


2.622 






di3 


3.068 


3.117 


3.050 


3.042 


3.039 


14 




d23 


2.879 


2.889 


2.853 


2.853 


2.853 


isomer 


I 


di2 


3.114 


3.144 


3.105 


3.105 


3.104 


14 




d 2 3 


2.888 


2.905 


2.868 


2.860 


2.854 


isomer 


II 


di2 


3.099 


3.155 


3.101 


3.090 


3.095 


14 




di4 


2.713 


2.735 


2.695 


2.680 


2.672 


isomer 


III 


d 2 3 


2.851 


2.857 


2.842 


2.823 


2.825 






di 5 


3.035 


3.053 


3.012 


3.012 


3.004 






d34 


3.353 


3.485 


3.360 


3.345 


3.352 


Li 6 + 




di2 


2.772 


2.800 


2.754 


2.734 


2.729 


isomer 


IV 


di4 


3.207 


3.250 


3.196 


3.186 


3.184 
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In view of these results, the intermediate basis set 3 
(6-31 1G), which yields results close to the ones obtained 
with the larger sets is chosen to perform the electronic 
structure calculations required in the MC simulations. 
The relative energies with respect to the most stable iso- 
mer, isomer IV, for isomers I, II, and III, are 17.58, 15.16, 
and 9.67 kJ/mol, respectively, for calculations using ba- 
sis set 3, compared with 19.38, 17.08, and 11.01 kJ/mol, 
respectively, using the large set 5. 

B. Electronic properties: dipole moments, Mulliken 
population charges and eigenvalues 

The Li4 cluster is non-polar (vanishing dipole moment) 
for symmetry reasons. Lis + is a charged system, so its 
dipole moment depends on the choice of the origin. How- 
ever, it is customary to evaluate the dipole moment using 
the center of charge as the origin, and in that case it pro- 
vides a useful indicator of the asymmetry of the charge 
distribution, and could also be experimentally relevant. 
Isomers I, and II are non polar, isomer IV is only slightly 
polar, but isomer III is considerably polar. The dipole 
moments of isomers III and IV, computed using basis set 
3, are 1.163 D »d 0.017 D, respectively. Mulliken pop- 
ulation chargesEB are also useful indicators of the charge 
distribution. Results obtained with basis set 3 for L14 and 
the four isomers of Lis + ar e shown in Table II. Significant 
differences are observed between different isomers, even 
between isomers I and II which are very similar both, 
geometrically and energetically. Both quantities provide 
useful indicators of isomerization during the MC simula- 
tions in the Lis" 1 " case. 

In addition, we present the two highest occupied and 
two lowest unoccupied molecular orbital energies for Li4 
and the four isomers of Li 5 + in Table II. The HOMO- 
LUMO gap is quite large in all cases, with values around 
1 eV. Therefore, it is unlikely that either thermal or quan- 
tum fluctuations will contribute to its closure. 

TABLE II. Mulliken populations and orbital energies of Li4 
and Li^~ computed with basis set 3. (in rA and eV). Atoms 
are labeled as in Figure 2. The two lowest energy unocuppied 
and two highest energy ocuppied orbital energies are given. 





Li4 


14 (!) 


Li 5 + (II) 


Li+ (III) 


14 (IV) 


qi 


0.2067 


0.5403 


0.0419 


0.2250 


0.3959 




0.2067 


0.1149 


0.2395 


0.1808 


0.3603 


q.'i 


-0.2067 


0.1149 


0.2395 


0.1808 


0.3603 


(14 


-0.2067 


0.1149 


0.2395 


0.2250 


-0.0583 


qs 




0.1149 


0.2395 


0.1878 


-0.0583 


ejv-i 


-3.9334 


-7.2535 


-7.2355 


-7.8535 


-8.1536 




-2.8986 


-6.6867 


-6.6703 


-6.5560 


-6.6603 


ejv+l 


-2.0169 


-5.3772 


-5.0774 


-5.4975 


-5.7065 


ejv+2 


-1.6806 


-5.0904 


-5.0774 


-5.4562 


-5.7002 



V. RESULTS OF THE PIMC-DFT SIMULATIONS 

A. Sampling strategy and convergence of the path 
integral with the degree of discretisation 

We used a simple Metropolis algorithm for the PIMC 
simulations with each trial move consisting of an attempt 
to move all the normal modes associated with all the par- 
ticles. Two different step-sizes were used: S c and S s . The 
maximum displacement of the center-of-mass was set by 
5 C , and that of the normal modes of order k - associated 
with a length scale Uk - by tik x S s . We have analyzed 
the possibility of introducing an additional convergence 
parameter k* , such that modes with k < k* or k > P—k* 
are moved with a relatively small step size, 5 S x <7fc, while 
those associated with small length scale fluctuations are 
moved by amounts proportional to 81 x Ok {81 > S s ). 
However, it turned out that, in this particular 
single step size 5 S for all values of k was efficient enough. 
This is often not the case when a large number of Trotter 
slices is used. The various parameters were adjusted to 
keep the overall acceptance ratio around 50% though oc- 
casional runs were used with acceptance ratios between 
40% and 60%. The same displacement parameter for the 
center-of-mass 5 C was used for both classical and quan- 
tum simulations. 

Simulations on the Li 2 dimer using a classical potential 
fitted to first-principles calculations were used to check 
the convergence of various properties with the degree of 
discretisation. 

Table III gives the PIMC results for the expectation 
value of the potential and kinetic energies using the prim- 
itive action and the fourth-order corrected form of the 
action. It can be observed that convergence to within 
the statistical error bars occurs with a Trotter number of 
just 4 when using the fourth-order correction, in contrast 
to 16 when using the primitive action. Errors in (V) arc 
an order of magnitude less than those in (K) . 



TABLE III. Convergence data for Li 2 at 100 K. Number 
of MC configurations is 5.12 million with acceptance ratios 
between 0.4 and 0.7. Error bars are given in brackets. First 
two columns correspond to results using the primitive approx- 
imation and the last two to results obtained using the fourth 
order correction to the effective action. Energies are expressed 
in K. 



m 


(V) 


(K) 




(V) 


(K) 


1 


-9521.78 (0.12) 


300.0 (0 


0) 


-9485.82 (0.14) 


335.7 (0.1) 


2 


-9491.74 (0.15) 


330.3 (0 


5) 


-9454.41 (0.15) 


366.2 (0.7) 


4 


-9465.76 (0.26) 


355.4 (1 


0) 


-9446.69 (0.31) 


372.8 (0.9) 


8 


-9453.39 (0.45) 


366.5 (1 


5) 


-9447.10 (0.32) 


371.8 (2.4) 


16 


-9449.79 (0.40) 


373.2 (7 


6) 


-9447.50 (0.44) 


371.1 (7.0) 
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This is typical of the relatively facile convergence of 
expectation values of operators diagonal in the coordi- 
nates as opposed to those than must be evaluated using 
thermodynamic estimators. In our experience, structural 
quantities such as pair distribution functions converge 
even faster than the potential energy. The increase in 
the error bars - at constant number of MC configura- 
tions - as the Trotter number is increased is also typical 
of PIMC simulations. Based on our tests with L12 we 
used a Trotter number of 4 at 100K and 8 at 50K for the 
larger clusters. 

We have performed a series of classical and quantum 
Monte Carlo simulations for Li4 and Li^ + clusters, with 
temperatures ranging from 50 K to 200 K. Each of 
them consisted of 10000-15000 Monte Carlo steps, pre- 
ceeded by around 1000 steps of thermalization. Wc 
stored atomic coordinates, energies, eigenvalues, Mul- 
liken population charges and dipole moment for each MC 
step, for later analysis. 

In the following we will concentrate only in struc- 
tural parameters, one-electron eigenvalues, and dipole 
moments, leaving aside other thermodynamical proper- 
ties which would need longer simulations to reduce the 
statistical error bars to useful values. 



B. Results for Li4 and Lig 

As mentioned in Section IV, L14 has a single, deep 
minimum at the 1 A singlet state, in the form of a pla- 
nar rhombus. Due to this fact the cluster is very rigid 
and thermal and/or quantum effects basically sample the 
PES around the minimum. The Li^" cluster constitutes 
a somewhat richer example due to the existence of sev- 
eral isomers. In particular, it is interesting to analyze the 
possibility of thermal activation and quantum tunneling 
between different regions of configuration space. In order 
to explore different regions of the PES, classical simula- 
tions were started from three different isomers, namely 
the ground state (isomer IV), and two higher energy con- 
figurations (isomers I and III). After a few thousand MC 
steps it was observed that the second and third simu- 
lations were attracted towards the ground state basin, 
showing that our MC strategy is quite effective in equi- 
librating and exploring configuration space, possibly be- 
cause interconversion barriers were low. Based on this, 
quantum simulations were started from the ground state 
(isomer IV) and, in order to facilitate the detection of 
possible tunneling behavior, also from the first excited 
isomer (isomer III) . Again, the isomerization towards the 
ground state was observed in this latter; the polymer 
moved as a whole, without showing any signature of tun- 
neling. Let us mention that none of the higher-energy 
isomers appeared again during the simulation, although 
structures slightly reminiscent of isomer III (the closest 
in energy to the ground state) were observed. In other 
words, Li^ appears to be unable to sample metastable re- 



gions of configurational space out from the ground state. 

Figure 4(a) shows the radial distribution function g(r) 
of Li4 in the classical case, and for different temperatures. 
It can be observed that the peaks are approximately cen- 
tered at the optimized zero-temperature distances. Tem- 
perature effects consist basically of broadening the peaks; 
the first of them, corresponding to the first neighbour 
shell, almost disappears above 200 K . In Figure 4(b) 
we show the effects of the quantum nature of the nu- 
clei by comparing simulations performed at 50 K using 
the classical and quantum schemes. It can be observed 
that quantum effects generate a pronounced broadening 
of the peaks, thus demonstrating the importance of their 
inclusion. 

In Figure 5 we show the radial distribution functions 
g(r) for Lig in the classical (a) and quantum (b) cases. 
The main panels contain the distribution averaged over 
all the 5 particles. Two groups of atoms can be identi- 
fied: a first one composed of three atoms more strongly 
bound, which form the central triangle of the bipyrami- 
dal structure IV, and a second one composed by the two 
external atoms. The upper inset shows the partial g(r) 
corresponding to these two groups in the classical case at 
100 K, and the lower inset in the quantum case at 50 K. 
No qualitative difference with Li4 can be observed. 

The trends discussed above also hold for the electronic 
properties, i.e. the distribution of one-particle eigenval- 
ues. As can be observed in Figure 6(a), the HOMO and 
LUMO eigenvalue distributions for Lig" exhibit signifi- 
cant broadenings upon temperature increase. It is to be 
remarked that the widths are different for different eigen- 
values, a fact that could be reflected in the temperature 
dependence of the optical photoabsorption spectrumEH. 
Similar broadenings can be observed in the lower panel, 
corresponding to the quantum and classical simulations 
for Li^" performed at 50 K. As advanced above, the 
minimum HOMO-LUMO distance never becomes smaller 
than about 0.5 eV, so that quantum effects cannot pro- 
mote an eigenvalue crossing which would result in a ma- 
jor modification of the electronic properties. 

Another important quantity is the electric dipole mo- 
ment, because of its experimental relevance. For Li4 the 
dipole moment vanishes at zero temperature due to sym- 
metry considerations. However, at finite temperature the 
cluster samples regions of the PES characterized by a fi- 
nite dipole moment, such that the mean value is non- 
vanishing. A more pronounced effect of the same type 
can be observed in the quantum simulations. The aver- 
ages computed using the classical and quantum schemes 
at 50 K are (0.17 ± 0.07) and (0.29 ± 0.15) D, for Li 4 and 
(0.14 ± 0.07) and (0.25 ± 0.12) D for Li+, respectively. 
The quantum dipole moment averages and distributions 
obtained at 50 K are similar to those obtained classically 
at about 100 K. 

It is interesting to note that where structural proper- 
ties are concerned, the overall effect of considering quan- 
tum nuclei is qualitatively similar to the effect of in- 
creasing temperature in the classical simulations. For 
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the closed-shell Li clusters considered in this work the 
classical temperature equivalent to the quantum system 
at 50 K is around 150 K. However, the dipole moment 
and eigenvalue distributions obtained using the quantum 
mechanical scheme at 50 K are qualitatively similar to 
those obtained within the classical scheme at about 100 
K . This points out an important difference between the 
two types of correlations involved in the statistical me- 
chanics of quantum systems, namely coherent quantum 
fluctuations as opposed to incoherent thermal fluctua- 
tions. These appear to behave in different ways according 
to the physical properties under consideration. 

Further characterization of the wave packet behaviour 
of the nuclei is provided by the imaginary-time corre- 
lation function, R 2 (t - t') =< \r(t) - r(t')\ 2 >, where 
0<r = t-t'< 13%. The value at r = /3K/2 is par- 
ticularly important because it gives an estimation of the 
quantum derealization of the nuclei. The derealization 
length, or gyration radius, is about 0.15 rA for both 
Li4 and Li^" at 50 K. The imaginary time correlation 
function for Li4 is shown in Figure 7, averaged for atoms 
1 and 2 (lower curve) , and for 3 and 4 (upper curve). 
It can be noticed that the two atoms which are more 
strongly bound (1 and 2 in Figure 2) are less delocalized 
than the other two, as may be expected. 

VI. CONCLUSIONS AND OUTLOOK 

We have introduced a novel method for simulating 
the statistical mechanics of quantum nuclei interacting 
through first-principles potentials, i.e. that derive di- 
rectly from the electronic structure. The scheme pre- 
sented and discussed here combines a path integral de- 
scription of the nuclear variables with an adiabatic, 
ground state, density functional description of the elec- 
tronic degrees of freedom. In the present scheme we have 
choosen a specific (NLDF-Gaussian) formulation to solve 
the electronic structure problem, but it is important to 
stress that any other implementation is perfectly valid 
and compatible with the present scheme, e.g. ab initio 
quantum chemical approaches like Hartree-Fock or MP2, 
and/or different localized (LCAO, LMTO) or extended 
(pseudopotential or augmented PW) basis sets. More- 
over, the present scheme is extremely simple to interface 
with any electronic structure code, since the only input 
needed to compute the statistical weigths are the (self- 
consistently) converged energy and forces. MD sampling 
schemes are more involved in this respect. 

We have shown the adequacy and the performance of 
this methodology by simulating quantum nuclear effects 
in the clusters Li4 and Lig". The number of imaginary- 
time slices needed to achieve convergence to a relative 
error of 0.5% in the nuclear kinetic energy is 4 at a tem- 
perature of 100 K, and 8 at 50 K. The same level of 
accuracy is obtained only with 16 slices (at 100 K) if the 
primitive approximation to the action is used. This rep- 



resents a gain of a factor of 4, which is very important 
due to the high computational cost of these simulations. 
The level of gain depends, however, on the shape of the 
potential energy surface sampled by the nuclei. 

The results presented here for the above clusters show 
that, at temperatures below 50 K, quantum nuclear ef- 
fects are crucial to account for their structural and elec- 
tronic properties. Pair correlation functions are quite 
broadened with respect to the classical counterparts, to 
a level that similar distributions would be obtained for an 
effective temperature of about 150 K if only thermal, and 
not also quantum, fluctuations were considered. The re- 
sults at T — 50 K can be considered to be representative 
of the ground state, as quantum nuclear effects largely 
overcome thermal motion. Electronic properties like one- 
electron eigenvalues and dipole moments show the same 
type of broadened distributions, although the effective 
classical temperature appears to be slightly lower, around 
100 K, instead of 150 K. This is to emphasize that quan- 
tum effects cannot be readily mimicked by adding extra 
thermal fluctuations, because the correlations involved 
are of a completely different character: quantum motion 
is coherent while thermal motion is incoherent. 

The simulation technology presented here opens up the 
possibility of studying the role of light nuclei, especially 
protons, in biological and chemical systems. The ad- 
vantage of describing the electronic component using a 
localized basis set (as Gaussians), as opposed to a PW 
basis set, for isolated clusters and molecules as well as 
for gas-phase reactions, is obvious because the vacuum 
around is easily taken into account. In addition, very of- 
ten reactive condensed-phase systems and large biological 
molecules do not need a full quantum description because 
the relevant chemical processes occur in a circumscribed 
region of space. For example, enzymatic reactions re- 
quire a first-principles electronic description only in the 
vicinity of the active site, while the rest of the system 
can be treated by means of classical force fields. There- 
fore, hybrid schemes that combine quantum and classical 
mechanical descriptions in different spatial regionsc3 are 
appropriate for these situations and, in conjunction with 
the present methodology, constitute a general computa- 
tional approach appropriate for studying the effects of 
nuclear derealization in the above situations. In fact, an 
insight into the problem of quantum hydrogen-bonding 
in water has already appeared in the literature^, thus 
signalling the time for a new and exciting area of multi- 
disciplinary research. 
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FIG. 1. A schematic representation of a potential energy 
surface with two minima separated by a barrier. Dashed lines 
represent different possible values of the zero-point-energy 
(ZPE): (a) ZPE is larger than the barrier (resonant regime), 
(b) ZPE is between the bottom of the higher well and the 
top of the barrier (tunneling regime), (c) ZPE is below the 
bottom of the higher well (classical regime). 

FIG. 2. The single isomer of Li4 (left panel), and the four 
isomers of Lig (right panel). Isomer IV is the most stable, 
followed by isomer III. Isomers I and II are quite higher in 
energy, and are almost degenerate. 

FIG. 3. Dependence of the energetics of Li clusters on the 
size of the basis set (BS). BS are ordered in increasing size. 



FIG. 4. Pair correlation functions g(r) for L14. The upper 
panel shows the classical results for temperatures of 50 (solid 
line), 100 (dashed line), and 200 (dotted-dashed line) K. The 
lower panel shows the quantum (dashed line) and classical 
(solid line) results for T = 50 K. 

FIG. 5. Pair correlation functions g(r) for Lig . The up- 
per panel shows the classical results for temperatures of 50 
(solid line), and 100 (dashed line) K. The lower panel shows 
the quantum (dashed line) and classical (solid line) results for 
T = 50 K. The insets show the partial g(r) for two groups 
of atoms, one forming the central triangle and the other con- 
sisting of the two external atoms. 

FIG. 6. Distribution of HOMO and LUMO one-electron 
eigenvalues for LiJ. The upper panel shows the classical re- 
sults for temperatures of 50 (solid line), and 100 (dashed line) 
K. The lower panel shows the quantum (dashed line) and 
classical (solid line) results at T = 50 K. 

FIG. 7. Root-mean-square of the imaginary-time corre- 
lation function for LI4 at 50 K plotted versus k (with 
< k = f3h/P < f3h), averaged over atoms 1 and 2 (lower 
curve) , and over atoms 3 and 4 (upper curve) . P is the Trotter 
number, 8 in this case. 
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